Report Experiment_Preprocessing

MYC-deficiency impairs the development of effector/memory T lymphocytes

Authors :Mathis Nozais1*, Marie Loosveld1,2*, Saran Pankaew1, Clémence Grosjean1, Noémie Gentil1, Julie Quessada1, Cyrille Mionnet1, Delphine Potier1@ & Dominique Payet-Bornet1@

*These authors contributed equally: Mathis Nozais, Marie Loosveld; @ Corresponding authors: Dr Delphine Potier and Dr Dominique Payet-Bornet

Link to article : (TO come)

This code is made to be running on Seurat 301v2 Docker. Data and explanation about this code are available at : https://github.com/mathisnozais/MycPten

Any questions on this analysis, please contact Mathis Nozais (nozais@ciml.univ-mrs.fr) or Delphine Potier (potier@ciml.univ-mrs.fr)

[1] “You are starting the analysis from count matrix obtain with CellRanger and Cite-seq-count”

Loading the first experiment

Looking to HTO distribution accross sample :

#HTO distribution
par(mfrow=c(1,1))
par(las=2)
par(mar=c(5,15,3,3))
barplot(rowSums(hto1), main = "Sequenced HTO distribution", horiz=TRUE)

rowSums(hto1)

##   Spleen-MP    Spleen-M Spleen-ctrl    Spleen-P   Thymus-MP    Thymus-M 
##      688432      206663      265412      662165      196262      117402 
## Thymus-ctrl    Thymus-P 
##      156856      562594

Demultiplexing HTO

Demultiplex cells based on HTO enrichment

Demultiplexing results

Cells classification

datatable(as.matrix(table(hashtag1@meta.data$MULTI_ID)),colnames = "Number of cells")

Violinplot (features)

VlnPlot(hashtag1,features = c("nFeature_RNA"),pt.size = 0.1, log = TRUE,  group.by = "MULTI_ID")

Violinplot (HTO counts)

VlnPlot(hashtag1,features = c("nCount_HTO"),pt.size = 0.1, log = TRUE,  group.by = "MULTI_ID")

tSNEs based on HTO

# Calculate a distance matrix using HTO
hto.dist.mtx <- as.matrix(dist(t(GetAssayData(object = hashtag1, assay = "HTO"))))

# Calculate tSNE embeddings with a distance matrix
hashtag1 <- RunTSNE(hashtag1, distance.matrix = hto.dist.mtx, perplexity = 100)

HTO margin

Tsne<-data.frame(
  tSNE_1 = hashtag1@reductions$tsne@cell.embeddings[,1],
  tSNE_2= hashtag1@reductions$tsne@cell.embeddings[,2],
  gene= hashtag1@meta.data$HTO_margin
)

HTO= hashtag1@meta.data$MULTI_ID
Max=max(hashtag1@meta.data$HTO_margin)
Min=min(hashtag1@meta.data$HTO_margin)
ggplot(Tsne,aes(x=tSNE_1,y=tSNE_2))+geom_point(aes(color=gene,shape=HTO))+
           scale_colour_gradient2(low = "blue",mid="orange",high="red",name="HTO margin",midpoint=(Max+Min)/2)+scale_shape_manual(values = c(15,16,17,18,19,20,21,22,23,24,25))

Ridge plots

Visualize enrichment for selected HTOs with ridge plots

RidgePlot(hashtag1, assay = "HTO", features = rownames(GetAssayData(hashtag1,assay = "HTO"))[1:2],ncol = 2, group.by = "MULTI_ID")

RidgePlot(hashtag1, assay = "HTO", features = rownames(GetAssayData(hashtag1,assay = "HTO"))[3:4],ncol = 2, group.by = "MULTI_ID")

RidgePlot(hashtag1, assay = "HTO", features = rownames(GetAssayData(hashtag1,assay = "HTO"))[5:6],ncol = 2, group.by = "MULTI_ID")

RidgePlot(hashtag1, assay = "HTO", features = rownames(GetAssayData(hashtag1,assay = "HTO"))[7:8],ncol = 2, group.by = "MULTI_ID")

Sample Information

The analysis will be run on the sample 1.

During the sample loading, we filter cells that do not pass the following filters.
Here are the description of those parameters in the Seurat CreateSeuratObject function:

• min.genes: Include cells where at least 200 genes are detected
• min.cells: Include genes with detected expression in at least 3 cells

a <- length(colnames(hashtag1@assays$RNA@data))
print(paste("After those filters, the remaining cell number is", a), quote = FALSE) 

## [1] After those filters, the remaining cell number is 4187

#add Exp1 cell identity
HTO_cr1 <- row.names(subset(hashtag1@meta.data, MULTI_ID == "Spleen-ctrl" ))
HTO_ct1 <- row.names(subset(hashtag1@meta.data, MULTI_ID == "Thymus-ctrl" ))
HTO_mr1 <- row.names(subset(hashtag1@meta.data, MULTI_ID == "Spleen-M" ))
HTO_mt1 <- row.names(subset(hashtag1@meta.data, MULTI_ID == "Thymus-M" ))
HTO_pr1 <- row.names(subset(hashtag1@meta.data, MULTI_ID == "Spleen-P" ))
HTO_pt1 <- row.names(subset(hashtag1@meta.data, MULTI_ID == "Thymus-P" ))
HTO_pmr1 <- row.names(subset(hashtag1@meta.data, MULTI_ID == "Spleen-MP" ))
HTO_pmt1 <- row.names(subset(hashtag1@meta.data, MULTI_ID == "Thymus-MP" ))
HTO_d1 <- row.names(subset(hashtag1@meta.data, MULTI_ID == "Doublet" ))
HTO_n1 <- row.names(subset(hashtag1@meta.data, MULTI_ID == "Negative" ))

HTO_thymus1 = c(HTO_ct1,HTO_mt1,HTO_pt1,HTO_pmt1)
HTO_spleen1 = c(HTO_cr1,HTO_mr1,HTO_pr1,HTO_pmr1)
HTO_identified1 = c(HTO_thymus1, HTO_spleen1)


# Create a Seurat object without doublet and unassigned cells (remove "negative", "doublet" & "nothing"))
clean.subset1 <- subset(x = hashtag1, cells = HTO_identified1)
VlnPlot(clean.subset1,features = c("nFeature_RNA", "nCount_RNA"),pt.size = 0.1, log = TRUE,  group.by = "MULTI_ID")

a <- length(colnames(clean.subset1@assays$RNA@data))
print(paste("After removing doublets and negative cells, the remaining cell number is", a), quote = FALSE) 

## [1] After removing doublets and negative cells, the remaining cell number is 3811

Mitochondrial percentage versus nFeatures

df<-data.frame(multi.id=Seurat1@misc$old_meta_data$MULTI_ID,percent.mito=Seurat1@misc$old_meta_data$percent.mito,nFeature_RNA=Seurat1@misc$old_meta_data$nFeature_RNA)
ggplotly(ggplot(df,aes(x=nFeature_RNA,y=percent.mito,color=multi.id))+geom_point())

UMAP

ggplotly(DimPlot(Seurat1, reduction = "umap", group.by = "MULTI_ID", do.label = TRUE, pt.size = 1))

Loading the second experiment

par(mfrow=c(1,1))
par(las=2)
par(mar=c(5,15,3,3))
barplot(rowSums(hto2), main = "sequenced HTO distribution", horiz=TRUE)

rowSums(hto2)

##   Spleen-MP    Spleen-M Spleen-ctrl    Spleen-P   Thymus-MP    Thymus-M 
##     3569690     1668124      929963     3251282      484276      749822 
## Thymus-ctrl    Thymus-P 
##      310711     1203239

Demultiplexing results

Cells classification

datatable(as.matrix(table(hashtag2@meta.data$MULTI_ID)),colnames = "Number of cells")

Violinplot (features)

VlnPlot(hashtag2,features = c("nFeature_RNA"),pt.size = 0.1, log = TRUE,  group.by = "MULTI_ID")

Violin plots (HTO counts)

VlnPlot(hashtag2,features = c("nCount_HTO"),pt.size = 0.1, log = TRUE,  group.by = "MULTI_ID")

tSNEs based on HTO

# Calculate a distance matrix using HTO
hto.dist.mtx <- as.matrix(dist(t(GetAssayData(object = hashtag2, assay = "HTO"))))

# Calculate tSNE embeddings with a distance matrix
hashtag2 <- RunTSNE(hashtag2, distance.matrix = hto.dist.mtx, perplexity = 100)

DimPlot(hashtag2, group.by = "MULTI_ID")

HTO margin

Tsne<-data.frame(
  tSNE_1 = hashtag2@reductions$tsne@cell.embeddings[,1],
  tSNE_2= hashtag2@reductions$tsne@cell.embeddings[,2],
  gene= hashtag2@meta.data$HTO_margin
)

HTO= hashtag2@meta.data$MULTI_ID
Max=max(hashtag2@meta.data$HTO_margin)
Min=min(hashtag2@meta.data$HTO_margin)

ggplot(Tsne,aes(x=tSNE_1,y=tSNE_2))+geom_point(aes(color=gene,shape=HTO))+
           scale_colour_gradient2(low = "blue",mid="orange",high="red",name="HTO margin",midpoint=(Max+Min)/2)+scale_shape_manual(values = c(15,16,17,18,19,20,21,22,23,24,25))

Ridge plots

Visualize enrichment for selected HTOs with ridge plots

RidgePlot(hashtag2, assay = "HTO", features = rownames(GetAssayData(hashtag1,assay = "HTO"))[1:2],ncol = 2, group.by = "MULTI_ID")

RidgePlot(hashtag2, assay = "HTO", features = rownames(GetAssayData(hashtag1,assay = "HTO"))[3:4],ncol = 2, group.by = "MULTI_ID")

RidgePlot(hashtag2, assay = "HTO", features = rownames(GetAssayData(hashtag1,assay = "HTO"))[5:6],ncol = 2, group.by = "MULTI_ID")

RidgePlot(hashtag2, assay = "HTO", features = rownames(GetAssayData(hashtag1,assay = "HTO"))[7:8],ncol = 2, group.by = "MULTI_ID")

Sample Information

The analysis will be run on the sample 2.

During the sample loading, we filter cells that do not pass the following filters.

Used parameters in the Seurat CreateSeuratObject function: * min.genes: 3 . Include cells where at least 3 genes are detected * min.cells: 200 . Include genes with detected expression in at least 200 cells

a <-length(colnames(hashtag2@assays$RNA@data))
print(paste("After those filters, the remaining cell number is", a), quote = FALSE) 

## [1] After those filters, the remaining cell number is 12990

#add Exp2 cell identity
HTO_cr2 <- row.names(subset(hashtag2@meta.data, MULTI_ID == "Spleen-ctrl" ))
HTO_ct2 <- row.names(subset(hashtag2@meta.data, MULTI_ID== "Thymus-ctrl" ))
HTO_mr2 <- row.names(subset(hashtag2@meta.data, MULTI_ID== "Spleen-M" ))
HTO_mt2 <- row.names(subset(hashtag2@meta.data, MULTI_ID== "Thymus-M" ))
HTO_pr2 <- row.names(subset(hashtag2@meta.data, MULTI_ID== "Spleen-P" ))
HTO_pt2 <- row.names(subset(hashtag2@meta.data, MULTI_ID== "Thymus-P" ))
HTO_pmr2 <- row.names(subset(hashtag2@meta.data, MULTI_ID== "Spleen-MP" ))
HTO_pmt2 <- row.names(subset(hashtag2@meta.data, MULTI_ID== "Thymus-MP" ))
HTO_d2 <- row.names(subset(hashtag2@meta.data, MULTI_ID== "Doublet" ))
HTO_n2 <- row.names(subset(hashtag2@meta.data, MULTI_ID== "Negative" ))

HTO_thymus2 = c(HTO_ct2,HTO_mt2,HTO_pt2,HTO_pmt2)
HTO_spleen2 = c(HTO_cr2,HTO_mr2,HTO_pr2,HTO_pmr2)
HTO_identified2 = c(HTO_thymus2, HTO_spleen2)

# Create a Seurat object without doublet and unassigned cells (remove "negative", "doublet" & "nothing"))
clean.subset2 <- subset(x = hashtag2, cells = HTO_identified2)
VlnPlot(clean.subset2,features = c("nFeature_RNA", "nCount_RNA"),pt.size = 0.1, log = TRUE,  group.by = "MULTI_ID")

a <-length(colnames(clean.subset2@assays$RNA@data))
print(paste("After removing doublets and negative cells, the remaining cell number is", a), quote = FALSE) 

## [1] After removing doublets and negative cells, the remaining cell number is 10309

Adding ADT

# Load in the UMI matrix
umi <- GetAssayData(object = clean.subset2, slot = "counts")

# Load in the ADT count matrix
raw.adt <- Read10X(PATH_ADT_DATA2, gene.column = 1)
adt <- raw.adt[c(1:6),]

rownames(adt) <- c("CD4","CD5","CD8","CD25","CD44","CD69")

#create an empty matrix containing NAs
Cell.list <- colnames(GetAssayData(object = clean.subset2[["RNA"]], slot = "data" ) )
ADT.list <- c(unique(rownames(adt)))
mat.adt <- matrix(nrow = length(ADT.list), ncol = length(Cell.list))
rownames(mat.adt) = ADT.list
colnames(mat.adt) = Cell.list

# Get cell barcodes detected by both RNA and ADT
joint_bcs <- intersect(colnames(umi),colnames(adt))
adt <- as.matrix(adt[,joint_bcs])

# Fill the empty matrix with values when existing
mat.adt[,joint_bcs]<-adt[,joint_bcs]

# Add ADT data as a new assay independent from RNA
clean.subset2[["ADT"]] <- CreateAssayObject(counts = mat.adt[,colnames(clean.subset2)])

# Normalize ADT data, here we use centered log-ratio (CLR) transformation
clean.subset2 <- NormalizeData(clean.subset2, assay = "ADT", normalization.method = "CLR")

#Scale
clean.subset2 <- ScaleData(clean.subset2, assay = "ADT")

ADT list :

print (rownames(adt))

[1] “CD4” “CD5” “CD8” “CD25” “CD44” “CD69”

Mitochondrial percentage versus nFeatures

df<-data.frame(multi.id=Seurat2@misc$old_meta_data$MULTI_ID,percent.mito=Seurat2@misc$old_meta_data$percent.mito,nFeature_RNA=Seurat2@misc$old_meta_data$nFeature_RNA)
ggplotly(ggplot(df,aes(x=nFeature_RNA,y=percent.mito,color=multi.id))+geom_point())

UMAP:

ggplotly(DimPlot(Seurat2, reduction = "umap", group.by = "MULTI_ID", do.label = TRUE, pt.size = 1))

Merging our two experiments

Load separate R object

You can load objects done with the code above. Or our object ?? (link )

Integrating the 2 seurat objects with seurat integration (cca)

a <- length(gene1)
b <- length(gene2)
c <- length(common_genes)

print(paste("We identified ",a," expressed in sample1 and ",b," expressed in sample2.",c,"are in common in this two set and will the integrated in the merged and corrected object."), quote = FALSE) 

## [1] We identified  13750  expressed in sample1 and  15815  expressed in sample2. 13631 are in common in this two set and will the integrated in the merged and corrected object.

UMAP:

Analysis part

Sample Information

The analysis will be run on the replicate 1 and 2.

During the sample loading, we filter cells that do not pass the following filters. We also filter cells that are detected as human/mouse multiplet using their barcodes.
Here are the description of those parameters in the Seurat CreateSeuratObject function:

• min.genes: Include cells where at least this many genes are detected
• min.cells: Include genes with detected expression in at least this many cells

a <- length(colnames(exp1.2.integrated))

print(paste("After those filters, and merging MYC_PTEN_01 and MYC_PTEN02 the remaining cell number is",a), quote = FALSE) 

## [1] After those filters, and merging MYC_PTEN_01 and MYC_PTEN02 the remaining cell number is 13912

UMAP and clustering parameter

Merge checking

HTO

Orig.idents

T-cell selection

According to T-cell markers we will exclude Cd3d low clusters: 13 (Bcells), 11, 14, 17 (monocytes/macrophages). According to T-cell markers we will exclude Cd3d/Cd3e low clusters: 11 (Bcells), 13, 17, 18, 19 (monocytes/macrophages), 16, 14 (ILC/NK).

Known RNA B and T markers

re clustering

T-cell umaps

HTO

clustering

END OF PREPROCESSING

We obtain the final object with clustering to start the analysis

Session Info

sessionInfo()

## R version 3.5.3 (2019-03-11)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 16.04.5 LTS
## 
## Matrix products: default
## BLAS/LAPACK: /usr/lib/libopenblasp-r0.2.18.so
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=C             
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] knitr_1.23         RColorBrewer_1.1-2 magrittr_1.5      
##  [4] dplyr_0.8.1        DT_0.6.2           gridExtra_2.3     
##  [7] kableExtra_1.1.0   plotly_4.9.0       ggplot2_3.1.1     
## [10] Seurat_3.0.1      
## 
## loaded via a namespace (and not attached):
##   [1] nlme_3.1-140        tsne_0.1-3          bitops_1.0-6       
##   [4] webshot_0.5.1       httr_1.4.0          sctransform_0.2.0  
##   [7] tools_3.5.3         R6_2.4.0            irlba_2.3.3        
##  [10] KernSmooth_2.23-15  lazyeval_0.2.2      colorspace_1.4-1   
##  [13] npsurv_0.4-0        withr_2.1.2         tidyselect_0.2.5   
##  [16] compiler_3.5.3      rvest_0.3.4         xml2_1.2.0         
##  [19] labeling_0.3        caTools_1.17.1.2    scales_1.0.0       
##  [22] lmtest_0.9-37       readr_1.3.1         ggridges_0.5.1     
##  [25] pbapply_1.4-0       stringr_1.4.0       digest_0.6.19      
##  [28] rmarkdown_1.12      R.utils_2.8.0       base64enc_0.1-3    
##  [31] pkgconfig_2.0.2     htmltools_0.3.6     bibtex_0.4.2       
##  [34] highr_0.8           htmlwidgets_1.3     rlang_0.3.4        
##  [37] rstudioapi_0.10     shiny_1.3.2         zoo_1.8-5          
##  [40] jsonlite_1.6        crosstalk_1.0.0     ica_1.0-2          
##  [43] gtools_3.8.1        R.oo_1.22.0         Matrix_1.2-17      
##  [46] Rcpp_1.0.1          munsell_0.5.0       ape_5.3            
##  [49] reticulate_1.12     R.methodsS3_1.7.1   stringi_1.4.3      
##  [52] yaml_2.2.0          gbRd_0.4-11         MASS_7.3-51.4      
##  [55] gplots_3.0.1.1      Rtsne_0.15          plyr_1.8.4         
##  [58] grid_3.5.3          promises_1.0.1      parallel_3.5.3     
##  [61] gdata_2.18.0        listenv_0.7.0       ggrepel_0.8.1      
##  [64] crayon_1.3.4        lattice_0.20-38     cowplot_0.9.4      
##  [67] splines_3.5.3       hms_0.4.2           SDMTools_1.1-221.1 
##  [70] pillar_1.4.0        igraph_1.2.4.1      future.apply_1.2.0 
##  [73] reshape2_1.4.3      codetools_0.2-16    glue_1.3.1         
##  [76] evaluate_0.13       lsei_1.2-0          metap_1.1          
##  [79] data.table_1.12.2   httpuv_1.5.1        png_0.1-7          
##  [82] Rdpack_0.11-0       gtable_0.3.0        RANN_2.6.1         
##  [85] purrr_0.3.2         tidyr_0.8.3         future_1.13.0      
##  [88] assertthat_0.2.1    xfun_0.7            mime_0.6           
##  [91] rsvd_1.0.0          xtable_1.8-4        later_0.8.0        
##  [94] survival_2.44-1.1   viridisLite_0.3.0   tibble_2.1.1       
##  [97] cluster_2.0.9       globals_0.12.4      fitdistrplus_1.0-14
## [100] ROCR_1.0-7